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ABSTRACT 

The presence of low mass, degenerate secondaries in millisecond pulsar binaries of- 
fers the opportunity to determine an age for the binary system independent of the 
rotational properties of the pulsar. To this end, we present here a detailed calcula- 
tion of the evolution of a grid of low mass (< 0.5 M©) helium core white dwarfs. 
We investigate the effects of different Hydrogen layer masses and provide results for 
well-known optical band-passes. We supplement the OPAL opacity calculations with 
our own calculations for low effective temperatures (T e g < 6000 K) and also provide 
fitting formulae for the gravity as a function of mass and effective temperature. In 
paper II we shall apply these results to individual cases. 
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1 INTRODUCTION 

Studies of binary pulsars (see Phinney & Kulkarni 1994 
for references) and close double degenerate systems (Marsh, 
Dhillon & Duck 1995) have discovered a number of low mass 
objects. The masses of these optically faint stars are esti- 
mated to be ~ 0.1 - 0.5 Mg. Thus, it is believed that these 
are helium core white dwarfs, which are not massive enough 
to ignite core helium burning and burn to carbon (Sweigart 
& Gross 1978; Mazzitelli 1989). Such stars arise in binaries 
because the progenitor is disrupted by Roche lobe overflow 
at some point in its natural evolution, losing its envelope and 
leaving behind a low-mass, degenerate helium star (Kippen- 
hahn, Kohl & Weigart 1967). 

Many millisecond pulsar binaries contain such dwarfs 
and, as such, their cooling ages offer an estimate of the 
age of the system, independent of the pulsar spin-down age 
and thus an interesting check of the traditional assumptions 
made about the ages of millisecond pulsars. 

Our aim in this paper is to present an extensive grid 
of cooling models covering the parameter space of low-mass 
helium core white dwarfs, much in the spirit of previous 
studies of the more ubiquitous carbon/oxygen core white 
dwarfs ( Wood 1992; DAntona & Mazzitelli 1990 and ref- 
erences therein). In paper II we apply these models to the 
observations of the low-mass binary pulsar population. 

In section ^ we briefly review the physical mechanisms 
which contribute to the white dwarf cooling process and 
discuss their relative importance. Section |§] describes our 
calculations of the low-temperature opacities necessary to 
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obtain accurate cooling sequences for the oldest systems. In 
section ^| we describe our numerical model and the tests of 
the code against other models from the literature. Finally, 
in section H we present our cooling sequences and describe 
the details of the cooling models. 



2 WHITE DWARF COOLING 

The basic qualitative picture of white dwarf cooling is a 
well-known one, going back to Mestel (1952). The star is 
supported by the pressure of a degenerate electron gas while 
the heat content is dominated by the thermal reservoir of 
non-degenerate ions, which can form a gas or a crystalline 
solid, depending on the Coulomb coupling parameter V — 
(Ze) 2 /akT, where a is the radius of the Wigner-Seitz sphere 
surrounding each ion. While the luminosity is dominated by 
the loss of thermal energy, young, hot white dwarfs may also 
have a contribution due to residual hydrogen burning at the 
base of the hydrogen envelope. 

The energy transport throughout most of the star is 
dominated by conduction due to the degenerate electrons, 
which keeps the core almost isothermal. The core is sur- 
rounded by a thin, non-degenerate envelope where the en- 
ergy transport is by radiative diffusion and, at late times, 
by convection (Bohm et al. 1977). This is the region of 
least efficient energy transport and thus it determines the 
rate at which the star cools. Hence, the chemical composi- 
tion of these outer layers is important. The high gravities 
of white dwarfs leads to gravitational settling (Schatzmann 
1958; Dupuis et al. 1992 and references therein) which re- 
sults in chemical separation in the white dwarf envelope and 
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hence the atmospheric opacity is dominated by the lightest 
species remaining on the surface, usually either hydrogen or 
helium. 

Since the cooling is determined by the radiative trans- 
port, the accuracy of the cooling curves will be dependant on 
the accuracy of the opacities we use. Over most of the cool- 
ing sequence we use the radiative opacities from the OPAL 
group (Rogers & Iglesias 1992), which are good down to 
T = 6000 K. For conductive opacities we use the results of 
Itoh et al. (1983) and Mitake, Ichimaru & Itoh (1984) for 
the region where the ions are gas or liquid and the electrons 
degenerate (helium cores don't reach crystallization temper- 
atures within a Hubble time). At temperatures and densi- 
ties not covered by the Itoh results, we use the older con- 
ductivities of Hubbard and Lampe (1969). This still leaves 
the region T < 6000 K unaccounted for. It is important be- 
cause many of the objects we will discuss in paper II have 
T e g ~ 4000 K. To address this issue, we have calculated 
Z = opacities for an arbitrary mix of hydrogen and he- 
lium using primarily the input physics of Lenzuni, Chernoff 
and Salpeter (1991), henceforth called LCS. In particular, we 
calculate the opacities from the collisionally induced absorp- 
tion (CIA) by molecular hydrogen. We describe this briefly 
below. 



3 LOW TEMPERATURE Z = OPACITIES 

We calculate the opacity of a gas containing H, He, H2, 
H _ ,He _ ,H + ,H^",H3",He + and electrons. We assume ideal gas 
and LTE. From the equation of state of Saumon, Chabrier 
and Van Horn (1995), we determine that the assumption of 
an ideal gas is good for p < 10 _2 gcm -3 . The partition func- 
tions and cross-sections are taken largely from LCS. For the 
calculation of collisionally induced molecular opacities we 
used the fits of LCS for the H2-H2 roto-vibrational tran- 
sitions (with corrections for typographical errors for which 
we thank Dr. Chernoff), after checking it with the original 
code (Borysow and Frommhold 1990), kindly supplied to 
us by Dr. A. Borysow. The roto-translational contribution 
was recalculated using another, more recent, code (Zhang & 
Borysow 1995) again supplied by Dr. Borysow. We also re- 
calculated the H2-He opacity following LCS although most 
of our atmospheres end up with pure H or He compositions. 

When we use the composition A"=0.72 and y=0.28, 
we find good agreement with LCS except at the high tem- 
perature, high density end T ~ 6000 — 7000 K and p ~ 
10 _3 gcm -3 , where we have a 10-20% error. We attribute 
this to our LTE treatment of the H~ ion. LCS note that 
the equilibrium abundance of this ion is affected by the ra- 
diation field at almost any temperature because of its low 
dissociation energy. Since the OPAL opacities reach down 
to 6000 K and our atmospheres are convective in this region 
anyway (see Figure ^|), we opt for simplicity and use our 
LTE results. 

For the case of a pure hydrogen atmosphere of moderate 
density p ~ 10 _5 gcm -3 , the opacity is dominated above ~ 
3000 K by the H~ ion and below by the collisionally induced 
opacity of H2. As noted by Bergeron, Saumon & Wesemael 
(1995), the Rosseland opacity goes through a minimum near 
this temperature. For a pure helium atmosphere He - and 



Rayleigh scattering provide most of what little opacity there 
is. 

We extend our calculations up to densities p ~ 1 g cm -3 . 
At higher densities, in lieu of an accurate calculation, we 
introduce abrupt pressure ionization and the opacity is de- 
termined by electron conduction at higher densities. For the 
case of a hydrogen atmosphere, this is unimportant because 
the atmosphere is sufficiently opaque that convection will 
dominate by the time these densities are reached. However, 
for the case of a helium atmosphere, this is the major source 
of uncertainty in our cooling times because the atmosphere 
is very optically thin and where we place the pressure ion- 
ization will determine the location of the photosphere. 

Figure [I] shows the contours of constant opacity in our 
density-temperature parameter space for pure hydrogen. It 
also shows the regions where different opacity tables were 
used. In Figure ^, as an example, we show the opacity as 
a function of temperature at p = 10~ 9 g cm -3 for a pure 
hydrogen atmosphere. There is excellent agreement between 
our calculation and the OPAL tables in the overlap region. 



4 THE COOLING CODE 

We have written a numerical cooling code using the Henyey 
method (Henyey et al. 1959,1964) to solve the stellar struc- 
ture equations. Our outer boundary conditions (obtained 
using a grey atmosphere calculation) are implemented using 
the Kippenhahn, Weigart & Hoffmeister (1967) method of 
triangles. 

Apart from the opacity tables described above, we have 
used the hydrogen and helium equation of state of Saumon 
et al. (1995), supplemented by a Thomas-Fermi model where 
necessary (at high densities). For the carbo n/ox ygen se- 
quences used in our comparisons (see section 4.1), we also 



used the carbon equation of state table of Fontaine, Gra- 
boske & Van Horn (1977). For completeness, we also include 
the effects of neutrino losses using the emissivities of Itoh & 
Kohyama (1983), Itoh et al (1984) and Munakata, Kohyama 
& Itoh (1985). To account for the residual nuclear burning, 
we incorporate the cross-sections of Fowler, Caughlan and 
Zimmerman (1975). We also use the mixing length theory 
of convection (we tested different parameterisations, ML1, 
ML2 and ML3 (see Bergeron, Wesemael & Fontaine 1992 
and references therein), and found no difference in their ef- 
fect on the cooling). Appendix |a] briefly describes the nu- 
merical details of the calculation. 



4.1 Code Tests 

The dearth of reliable helium white dwarf cooling sequences 
means that, in order to test our numerical method properly, 
we need to include carbon cores so as to test our models 
against those in the accepted literature. We have made com- 
parisons with three well known codes, adjusting our input 
physics to approximate the original input physics as closely 
as possible. The wide variation in input physics and parame- 
ters makes this a non-trivial proposition (see Winget & Van 
Horn 1987). Nevertheless we obtain satisfactory fits that give 
us confidence in the accuracy of our numerical scheme and 
also incidentally reinforces the assertion of Winget and Van 
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Horn that the wide variation in published white dwarf cool- 
ing models is a result of input physics and not numerical 
treatment. The comparison is shown in Figure ^. The three 
models are from D'Antona & Mazzitelli (1989), Koester & 
Schonberner (1986) and Wood (1992). Of particular interest 
is the comparison with the Wood cooling sequence because 
that is the most up-to-date, using Z = OPAL opacities. 
We don't expect agreement at late times because the cal- 
culations of the opacities for T < 6000 K will vary, but the 
agreement for ages less than a few billion years is excellent. 
In particular, the arrows indicate the point at which crys- 
tallisation begins, and we do not expect agreement at later 
times than this because our treatment of crystallisation is 
not sufficiently detailed for these purposes (Helium cores do 
not crystallise in a Hubble time and thus it is unimportant 
for this analysis). An updated version of these curves will 
presented in a future publication. 

The one useful test of our code for the case of a helium 
white dwarf is a comparison with the cooling sequence of 
Iben & Tutukov (1986), shown in Figure^. Their model was 
for a 0.3 Mq star which, after undergoing two hydrogen shell 
flashes, has a shell of 2.5 x 10 -4 Mq of pure hydrogen during 
the white dwarf phase. They caution that their opacities 
are uncertain at low luminosities, but it provides a useful 
comparison at least for logL/L© > —3. Their model also 
demonstrates another uncertainty for our higher luminosity 
models. The white dwarf stage really only begins after the 
end of the second hydrogen shell flash which means that our 
t = corresponds to t = 10 s years in this case. The effect 
of this uncertainty in the starting point for the model is 
important for logL/LQ > —3. 

A second uncertainty in the initial conditions is the 
value of the core temperature at the beginning of the evolu- 
tion. Tests of our code show that an uncertainty of 3 x 10 7 K 
corresponds to an age uncertainty ~ 10 8 years, although the 
exact numbers are model dependant. We solve this problem 
by using the results of evolutionary calculations (Soberman, 
pers. comm.) of the progenitors at a variety of binary sep- 
arations around a neutron star, using the stellar evolution 
code of Eggleton (1971) (see also Pols, Tout & Eggleton 
1995). This provides us with an estimate of the initial cen- 
tral temperature for a given remnant mass, which we use as 
a starting point for our calculations. For most of the sys- 
tems that we address in paper II, with ages > 1 Gyr, this 
uncertainty is unimportant. 

There are also uncertainties in the evolution associated 
with the uncertainties in the input physics. We investigated 
the effect of extrapolating the Itoh opacities to lower densi- 
ties instead of using the old Hubbard and Lampe opacities 
for those regions. This was found to have little impact on 
the cooling sequences, because the conductivity is most im- 
portant in the highest density regions, which is covered by 
the Itoh tables. A second uncertainty is the effect of extrap- 
olating the OPAL opacities to higher densities in the region 
T ~ 10 4 — 10 5 K and p ~ 10~ 3 - 1 gem" 3 . This is not terri- 
bly important because the atmosphere is convective in this 
region, and so the extrapolation will affect only the depth 
at which convection terminates. 

A major uncertainty of past calculations was for the 
low temperature (T < 6000 K) opacities. Our calculations 
in section ^| are designed to solve this problem, at least for 
the hydrogen atmospheres. The problem persists with the 



helium atmospheres, where the photospheric opacity is de- 
termined by the pressure ionization of helium, in a regime 
where accurate opacity calculations do not exist. For some 
models, we also have convective mixing of the elements. We 
assume that the convective zone is mixed to a uniform com- 
position determined by the relative mass fractions of hydro- 
gen and helium in the convective zone. When the convective 
zone retreats from the lower layers, we assume that the sep- 
aration is instantaneous, i.e., at all times, the composition in 
the convective zone is determined by the depth of the con- 
vective base (i.e., the mass of helium in the convection zone 
relative to the hydrogen envelope mass). The admixture of 
helium into the hydrogen atmosphere does not restrict the 
validity of our opacity calculations at low temperatures. Be- 
cause of the extremely low neutral helium opacities, 5% of 
hydrogen by mass is still sufficient to provide enough opacity 
so that the photospheric pressure lies well below the pressure 
ionization value. 



5 RESULTS 

In this section, we will describe in detail the cooling of a 
0.3 Mq helium core star with a thick (3 x 1O _4 M0) hydro- 
gen envelope, shown in Figures |E] to [?]. We start our models 
with an age of 10 s years (the time taken for the Iben & Tu- 
tukov models to reach the end of the last shell flash). We 
find that, even for these hydrogen masses and at early ages, 
the residual nuclear burning contribution to the luminosity 
is never more than ~ 1% (see Figure^). The neutrino lumi- 
nosity is never more than ~ 10% of the photon luminosity. 
Thus, over the first few xlO 6 years the star completes the 
contraction to the white dwarf configuration that it began 
when the last shell flash ended (see, e.g., Iben & Tutukov 
1986, Table 1). In these early stages, the helium core is only 
mildly degenerate (central degeneracies ~ 30) with r\ ~ 1 at 
the hydrogen-helium boundary, although the core is already 
approaching an isothermal state (0.99 of the stellar mass has 
T > 0.5 T c ). The hydrogen envelope is entirely radiative at 
this stage, so that this is the time at which the star most 
resembles the Mestel ideal (see Figure |^) . 

When r efi ~ 15000 K (~ 3 x 10 s years), a small con- 
vection zone appears near the surface. This zone remains 
relatively thin until Tea ~ 10 4 K ( ~ 10° years), when it 
starts to deepen as a consequence of the movement of the 
hydrogen ionization zone to greater depths (Figure . The 
recombination of hydrogen also leads to an increase in the 
photospheric density (since neutral hydrogen has a smaller 
opacity than ionized hydrogen) . This increase is only halted 
once molecular hydrogen begins to form (T e a ~ 5000 K, 
t ~ 4 x 10 9 years, see Figure H) . These photospheric changes 
also lead to the flattening of the cooling curves and the in- 
crease in the gravity seen in Figures [l^ and [IB], similar to 
effects in C/O models as described in D'Antona & Mazzitelli 
(1990). The deepening of the convective zone continues until 
the base reaches a depth where the conduction due to degen- 
erate electrons is more efficient than convection. This occurs 
after 3 x 10 9 years. At its deepest extent, the convection zone 
contains ~ 5 x 1O _5 M0 . This is only a couple of scale heights 
above the hydrogen/helium interface, so we might expect 
some small amount of photospheric helium contamination 
(since the convection extends to optical depths r ~ 0.3) be- 
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low T e g ~ 6000 K (see also Bohm et al 1977). However, the 
helium will be present only in trace amounts and will not af- 
fect the cooling. As the star cools further, it will revert again 
to a pure hydrogen atmosphere because of the formation of 
a radiative buffer zone near the surface below T e g ~ 3500 K. 
This occurs because there is a minimum in the atmospheric 
opacity as the primary opacity contribution changes from 
CIA of H2 to H~ absorption. If this minimum value is low 
enough, the region of the minimum will be radiative, but 
convection will continue to operate both above and below. 
This would allow the surface helium to diffuse out of the 
top convection zone, leaving the atmosphere once again in 
a pure hydrogen state. 

The features of the above sequence are common to most 
of the models we discuss here. The primary differences occur 
at the lower mass end of the models (M ~ 0.15 Mq). The 
lower masses lead to lower central densities and lower degen- 
eracy, so that the model radii are somewhat larger (hence 
with lower gravities) and convection extends deeper into the 
cooler models (because conduction is less efficient). Thus, 
even for the thickest hydrogen envelopes considered hereQ 
the 0.15 Mq models do dredge up some helium, becoming 
mixed atmosphere stars for a short while until the base of 
the convective zone retreats again. Figure |^ shows a dia- 
gram similar to Figure ^ but for a 0.15 Mq model. We see 
that the base of the convective zone penetrates to the he- 
lium layer, leading to an atmospheric helium abundance of 
~ 15% in this case. Hence, for about 1.5 Gyr, the star would 
exhibit the characteristics of a cool H/He atmosphere star, 
although, with a temperature of 4000-5000 K, it would be 
difficult to identify it as such. For the thick hydrogen layer 
models, only the lowest mass models are affected by atmo- 
spheric helium contamination. In the sequence of models 
with thin (~ 10 _6 Mq) hydrogen masses, dredge-up occurs 
for all masses. In Figure M we show the effective temperature 
ranges in which this occurs for all model masses. 

Figures and show the cooling sequences for six dif- 
ferent models spanning the range of representative masses. 
We show curves for 0.15, 0.25 and 0.45 Mq, with two dif- 
ferent hydrogen envelope masses, 3 x 10 _4 Mq (thick) and 
10 _6 Mq (thin). For completeness we also show the se- 
quences for pure helium models in Figure [l^. In principle, 
as found by DAntona & Mazzitelli (1989), one can obtain 
shorter cooling times for these models, but the opacities for 
T e g < 6000 K are rather uncertain, and thus we do not ex- 
pect the cooling curves to be accurate at temperatures below 
this value. 

In order to make comparison with observations easier, 
we have also calculated the black-body absolute magnitudes 
for these sequences, using the flux calibrations of Bessell 
(1979). The most common bands used are V and I. By com- 
parison with the detailed atmospheric calculations of Berg- 
eron, Saumon & Wesemael (1995), we find that this is an ex- 



t We have used the upper limit derived by Iben & Tutukov 
(1986), although it is likely to increase somewhat for smaller core 
masses. It is important to note at this point that many stellar 
evolution calculations such as those by Pylyser & Savonije (1988) 
or Alberts et al (1996) cannot be used to set this limit, because 
they make use of a stellar evolution code which cannot resolve 
the shell flashes described by Iben & Tutukov. 



cellent approximation for T e s > 4500 K. This can be under- 
stood by comparing the Rosseland & Planck mean opacities 
shown in Figure |^ which correspond to different frequency 
weightings of the true monochromatic opacity. When they 
are equal this indicates an approximately uniform contin- 
uum opacity and the black body approximation is valid. For 
Tcff < 4500 K, the CIA of H2 begins to dominate and this 
leads to significant deviations from the black body approxi- 
mation, as can be seen in Figure ^. Although the Bergeron et 
al results are calculated for higher gravities than we require, 
they represent an upper limit on the deviations from a black- 
body. This is because CIA is an intrinsically high density 
phenomenon (it requires collisional interactions) and pho- 
tospheric densities are lower in lower gravity atmospheres. 
This can be demonstrated by comparing the low tempera- 
ture results for the log g = 7.5 and log g = 8.5 atmospheres 
of Bergeron et al. We will discuss the uncertainties again in 
paper II when we come to apply these models. 

In Tables |l| ^ and ^| we present the model parameters 
(mass, age and gravity) as a function of the observable quan- 
tities Teff and Mi. Each entry contains three numbers, the 
mass in Mq, the age in Gyr and log g. Figures [[3] and |l4| 
show My and Mi as a function of age. 

For stars with hydrogen envelopes and effective temper- 
atures > 7000 K, spectroscopic determinations of effective 
temperature and gravity can provide a direct measurement 
of the white dwarf mass (Bergeron, Wesemael & Fontaine 
1991) , provided one has a relationship between radius and 
mass. We have calculated the mass-radius relations for both 
our thick and thin hydrogen layer models for the full range 
of helium core white dwarf masses. We fit this by a relation 
between gravity g = GM/R 2 and effective temperature T e ff 
for a given mass, namely 



g = F(Teff) [02 - a-jTeff} , 
with 

F(Teff) = 1 + 



aiT e f[ 



1 + 9 exp (4 x IO- 3 (T off - 5800K)) 1 



and 



01 = 9.91 x 10' 

a 2 = 1.69 x 10 
M 



a s = 946 



7 ( M_ 

Mq" 

, 1 M_ 

Mq 
0.17 



M 



for the thick hydrogen envelope and 

-2.18 

- / W \ 

ai = 3.79 x 10' 



7 I M 

Mb 



a 2 = 1.68 x 10 s 



M 

Mq" 



(1) 

(2) 

(3) 
(4) 
(5) 

(6) 

(7) 
(8) 



as = 650 

for the thin hydrogen envelope and M < 0.4 Mq. Figure |l5| 
compares the fits to the proper curves for the case of a thick 
envelope. 

We may compare the T e g = limit of ([j]) with vari- 
ous well-known T = mass-radius relations. Converting the 
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above into mass-radius relations and extrapolating to the 
T — limit we get 



i| = °-° 13 (l£ ) ( ™ CkH) ' 



(9) 



(10) 



This is very close to the often used mass-radius relationship 
of Paczynski (1967) for low mass degenerate dwarfs P/Rq = 
0.013 (M/M Q )~ 1/3 . In Figure |l| we compare our T — 
curves with the Paczynski relation and that of Hamada and 
Salpeter (1961). 



6 CONCLUSION 

We have presented a set of cooling sequences for low mass 
helium white dwarfs of different masses and with different 
masses of surface hydrogen. We provide blackbody absolute 
magnitudes and surface gravity - effective temperature re- 
lations as an aid to the analysis of future observations. In 
paper II we shall apply these models to the optical obser- 
vations of the companions to millisecond pulsars in order to 
derive cooling ages. 

The authors would like to thank Marten van Kerkwijk 
and Yanqin Wu for lengthy discussions about white dwarf 
physics and Glenn Soberman for helping with the initial con- 
ditions for our models. The generosity of messrs D. Saumon, 
G. Fontaine, I. Mazzitelli, F. Rogers and C. Iglesias in pro- 
viding their microphysical results is also appreciated. We 
would also like to thank the referee, Dr. D'Antona, for some 
insightful comments on the atmospheric physics. 



APPENDIX A: BRIEF NUMERICAL DETAILS 

It is our intention to make both the results and the cool- 
ing code publicly accessible. To obtain the most up-to-date 
version of the code email hanaen@cita.utoronto.ca. In this 
brief appendix we outline the version of the stellar structure 
equations that our program is designed to solve. 
The four stellar structure equations are 



dP GM(r) 
dr r 2 
AM „ 2 

— — = 47IT p, 

dr 

dL 2 

__ — 47T)- p e 

dr 



dT 
dP 



T. 



where V in equation (A.4) is determined by 



V = 



L(r) 



ISnac T 3 GM(r) 



(Al) 
(A2) 
(A3) 
(A4) 

(A5) 



for a radiative atmosphere or from a mixing length calcu- 
lation for a convective atmosphere. The function e contains 
both the nuclear energy generation terms and the thermal 
and gravitational loss terms, as well as the neutrino energy 



losses. We use the transformation of variables from (see Ta- 
ble k| for definition of symbols) 



x = r/10 11 cm, 

t = L/Lo, 

7T = InP/Po, 

6 = lnT/To, 

<t> = hip/po, 
L 
M 
Lppo 
T ~ M Po ' 
m = M/M = 1 - exp(-£), 

with the normalizations 



Po 
Po 

T 

and 
L 



2.123 x 10 14 dyncm~ 2 



0.159 gem 



M 
M, 







3.568 x 10 6 K | i^P^ 



Mp Lp 
M L 



Mo 

M 



1/4 



(A6) 
(A7) 
(A8) 
(A9) 
(A10) 

(All) 

(A12) 
(A13) 

(A14) 
(A15) 

(A16) 

(A17) 
(A18) 



100 L 

The set of equations that we actually solve numerically is 
d£ _ 

d£ x 2 ' 

d£ e / Q + r dTr _ dfli 

d£ - Wad l ad dr dr 



Wad 
€-v ~\~ ^nuc) j 



(A19) 

(A20) 
(A21) 

(A22) 



dC 1 L e I x i ' 

d£ v ' x 4 

(Q = - (flff ) p , Vad = (HS) S )- This s y stem of equations 
is solved on a grid equally spaced in £. The external bound- 
ary condition is determined by integrating inwards using an 
adaptive stepsize algorithm starting with a grey atmosphere 
solution to determine n and 9 at the outer grid point as a 
function of £ and x. The fitting procedure used was taken 
from the triangle method of Kippenhahn et al. (1967). 
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Table 1. Masses of helium Core white dwarfs (Mq) 

T e ff 
(K) 



10.0 


10.5 


11.0 


11.5 


12.0 


12.5 


0.27 


0.36 










0.25 


0.33 


0.45 








0.24 


0.30 


0.43 








0.22 


0.28 


0.40 








0.21 


0.26 


0.36 








0.19 


0.24 


0.32 


0.45 






0.18 


0.22 


0.29 


0.42 






0.17 


0.20 


0.26 


0.38 






0.16 


0.18 


0.23 


0.32 








0.16 


0.20 


0.27 


0.42 








0.18 


0.23 


0.34 








0.15 


0.19 


0.27 


0.43 








0.16 


0.21 


0.33 










0.19 


0.29 










0.17 


0.25 










0.15 


0.21 












0.18 



12.75 13.0 



Mi 
13.25 



13.5 13.75 14.0 14.25 14.50 14.75 15.0 



12000 

11500 

11000 

10500 

10000 

9500 

9000 

8500 

8000 

7500 

7000 

6500 

6000 

5800 

5600 

5400 

5200 

5000 

4800 

4600 

4400 

4200 

4000 

3800 

3600 

3400 



0.42 
0.37 
0.32 
0.27 
0.22 
0.19 
0.15 



0.44 
0.41 
0.35 
0.29 
0.24 
0.20 
0.16 



0.44 
0.39 
0.32 
0.26 
0.20 
0.16 



0.42 
0.34 
0.27 
0.21 
0.16 



0.44 
0.37 
0.28 
0.22 
0.16 



0.39 








0.29 


0.40 






0.22 


0.29 


0.40 




0.16 


0.21 


0.29 


0.40 






0.20 


0.27 








0.18 



0.37 
0.24 



The values in this table are for a thick hydrogen envelope (3 X 10 4 Mq). 



Table 2. Ages of helium core white dwarfs (Gyr) 



Teff 

(K) 



10.0 10.5 11.0 11.5 12.0 12.5 12.75 13.0 



Mi 
13.25 



13.5 13.75 14.0 14.25 14.5 14.75 15.0 



12000 

11500 

11000 

10500 

10000 

9500 

9000 

8500 

8000 

7500 

7000 

6500 

6000 

5800 

5600 

5400 

5200 

5000 

4800 

4600 

4400 

4200 

4000 

3800 

3600 

3400 



0.2 
0.2 
0.2 
0.2 
0.2 
0.1 
0.1 
0.1 
0.1 



0,1 
0,1 
0.4 
0.4 
0.4 
0.1 
0.3 
0.3 
0.3 
0.3 



0.7 
0.8 
0.8 
0.7 
0.7 
0.6 
0.6 
0.6 
0.5 
0.5 
0,1 



1.1 
1.2 
1.1 
1.1 
1.0 
0.9 
0.8 
0.7 



1.8 














1.6 














1.5 


2.6 












1.3 


2.3 


3.1 










1.2 


2.1 


2.9 


3.9 








1.1 


2.0 


2.6 


3.7 








1.0 


1.8 


2.4 


3.4 


4.9 








1.6 


2.2 


3.1 


4.6 










2.0 


2.8 


4.0 


5.9 








1.8 


2.5 


3.5 


5.1 


7.1 








2.2 


3.1 
2.7 


4.4 
3.7 
3.1 


6.1 

5.3 
4.3 
3. 1 



7.7 








6.1 


8.9 






4.8 


6.9 


10.0 




3.7 


5.3 


7.6 


10.9 






5.5 


7.9 








5.5 



11.3 
7.9 



The values in this table are for a thick hydrogen envelope (3 X 10 4 Mq). 
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Table 3. Gravities of helium core white dwarfs (log g) 



(K) 



10.0 


10.5 


11.0 


11.5 


7.01 


7.33 






6.95 


7.27 


7.60 




6.88 


7.19 


7.54 




6.81 


7.12 


7.47 




6.74 


7.04 


7.38 




6.66 


6.95 


7.28 


7.62 


6.57 


6.85 


7.18 


7.55 


6.50 


6.75 


7.07 


7.44 


6.39 


6.65 


6.95 


7.30 




6.53 


6.82 


7.16 






6.68 


7.00 






6.53 


6.83 








6.65 



12.0 12.5 12.75 13.0 



Mi 
13.25 



13.5 13.75 14.0 14.25 14.5 14.75 15.0 



12000 

11500 

11000 

10500 

10000 

9500 

9000 

8500 

8000 

7500 

7000 

6500 

6000 

5800 

5600 

5400 

5200 

5000 

4800 

4600 

4400 

4200 

4000 

3800 

3600 

3400 



7.54 
7.37 
7.18 
6.97 
6.87 
6.78 
6.67 



7.57 
7.35 
7.25 
7.13 
7.01 



7.56 
7.45 
7.34 
7.22 
7.07 
6.93 
6.77 



7.64 
7.55 
7.43 
7.29 
7.14 
6.98 
6.81 



7.62 
7.51 
7.36 
7.20 
7.03 
6.84 



7.58 
7.42 
7.25 
7.07 
6.86 



7.63 
7.48 
7.29 
7.08 
6.86 



7.52 








7.31 


7.55 






7.09 


7.32 


7.56 




6.84 


7.07 


7.30 


7.55 






7.02 


7.26 








6.94 



7.50 
7.18 



The values in this table are for a thick hydrogen envelope (3 X 10 4 Mq). 



Table 4. Symbols used in Appendix 
Physical quantity Symbol Transformed 



Pressure 


P 


77 


Temperature 


T 


e 


Luminosity 


L 


t 


Radius 


r 


X 


Mass 


m 




Time 


t 


r 


Density 


P 
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Figure 1. The Hydrogen Phase Diagram: The thin solid 
lines are contours of constant log re, where re is the Rosseland mean 
opacity in cm 2 g _1 . The contours have values 0, 2 and 4. The 
dotted lines are also contours of constant log re but with values 
-2, -4, etc., down to -12 (decreasing monotonically in all directions 
from the opacity peak near 30000 K and 10 — 3 g.cm~ 3 ). The heavy 
dashed lines delineate regions where different tables have been 
used to calculate the opacity. The Itoh opacities are valid for T < 
0.1 Tp (where Tp is the Fermi temperature), y < 0.1 [y measures 
the importance of the wave nature of the ions) and T < 171. There 
is also a lower bound on the density p > 100 gem -3 . Outside of 
this region we use the conductivities of Hubbard and Lampe. The 
radiative opacities in the region T ~ 10 4 — 10 5 K and p ~ 10 — 4 — 
lgcm -3 were obtained by extrapolating the OPAL opacities to 
higher densities. This extrapolation is not important because the 
atmosphere is convective at these temperatures and densities. The 
box in the lower left-hand corner is the region covered by our 
opacity calculations. We can see the opacity minimum near 3000 
K due to the change in the dominant opacity mechanism from 
H — absorption to H2 CIA. The strange behaviour in the upper 
left-hand corner is due to the extrapolation of the conductive 
opacities outside their range of validity. This is unimportant as 
no model we consider will approach this region. 



Figure 2. Matching Opacities: The solid points are Rosseland 
mean opacities calculated using our code. The open squares are 
the results of the OPAL calculation. The open stars are again 
our calculation but showing Planck mean opacities. The vertical 
dotted lines delineate the region 6000-7000 K which is where the 
two calculations overlap. Once again, the minimum in the opac- 
ity near 3000 K is due to the change in the dominant opacity 
contributor, from H2 CIA at lower temperatures to H — absorp- 
tion at higher temperatures. We see that below 5000 K there is a 
significant discrepancy between the Planck and Rosseland mean 
opacities. 



Figure 3. Code Comparisons: C/O models: The open 
squares denote the models we compare against. The filled squares 
are our own models. The arrows denote the point at which the 
core of the model begins to crystallise, although some of the dif- 
ference at low luminosities is due to the updated opacities. The 
left panel describes a 0.564 Mq oxygen core surrounded by a 
helium envelope of 2.5 X 10 — 3 Mq and a hydrogen envelope of 
3 X 1O~ 4 M . The metallicity is taken to be Z = 0. The centre 
panel describes a 0.546 Mq carbon core with a helium envelope 
of 0.022 Mq and hydrogen envelope 10 _4 Mq. The metallicity is 
Z = 0.02. The rightmost panel is a 0.6 Mq star, with a carbon 
core, mass fraction 10 — 2 of helium and 10 -4 of hydrogen. The 
metallicity is Z = 0. We had to adjust our conductive opacities 
to reproduce the above results. When prior authors used Hubbard 
& Lampe opacities in regions where we used Itoh opacities, we 
divided our opacities by a factor of 2 to compensate. 



Figure 4. Code Comparisons: He model: The solid line 
shows our model with no corrections for different starting points. 
The dashed line corresponds to the same model, but with the 
age incremented by 10 8 years, to compensate for the time spent 
in prior evolutionary stages. The open squares are the results of 
Iben & Tutukov (1986). The agreement is excellent until t ~ 2 
Gyr, by which point Iben & Tutukov caution that their opacities 
are uncertain. 



Figure 5. Cooling of a 0.3 Mq Model: We show here the cool- 
ing of a 0.3 Mq star with a hydrogen envelope of 3 X 10 _4 Mq. The 
solid line is the electromagnetic bolometric luminosity Lj, the 
dotted lines indicate the neutrino and nuclear luminosities respec- 
tively (the nuclear contribution is included in the bolometric lu- 
minosity) and the short and long dashed lines indicate the gravity 
(in units of lO'cms -2 ) and the central degeneracy (r/ c = Ep/kT) 
respectively. The effective temperatures corresponding to the var- 
ious ages for this model are shown on the top axis, ranging from 
20000 K at the left to 3000 K at the right. 



Figure 6. Evolution in the Phase Diagram: This T—P phase 
diagram shows three representative atmosphere profiles (labelled 
at the top by their age in Gyr) for the evolution of the same 
model shown in |i[ The heavy solid lines indicate the helium parts 
of the star and the thin solid lines indicate the hydrogen part. The 
dotted lines delineate the regions of 50-50 division between HI-HII 
and H2-HI respectively. The dashed line indicates the boundary of 
the convective region for this model (the other pair of dashed lines 
in the upper left-hand corner indicates crystallization boundary 
of helium). The thick solid line at the lower left indicates the 
location of the photosphere for this cooling sequence. The labelled 
dashes on each of the three curves indicate the points at which 
the degeneracy parameter rj = Ep jkT has that particular value. 



Figure 7. 0.3 Mq Convective Zone: Here we show the mass 
in the convective zone as a function of age (or effective tempera- 
ture). The shaded region is the convective zone, while the heavy 
solid line indicates the location of the photosphere. The dashed 
line indicates the hydrogen-helium interface in this set of models. 
We note the appearance of a radiative buffer zone at late times 
associated with the transition from H _ opacity to H2 opacity. 



Figure 8. 0.15 Mq Convective Zone: Once again the shaded 
region is the convective zone, and the heavy solid line denotes the 
position of the photosphere. We note that, for T c fj ~ 4000 — 5000 
K, the atmosphere will be contaminated with helium. 



Figure 9. Atmospheric Helium Contamination for Thin 
Hydrogen Envelopes: The shaded regions indicate those mod- 
els in which the convective zone extends into the helium layer and 
thus causes atmospheric helium contamination. We consider two 
representative cases. The area marked as H/He is characterised by 
a mass fraction of helium, Xh c > 0.1. The area marked as He/H 
is characterised by Ajj c > 0.8. This cutoff value may seem rather 
high, but we note that the dredge-up occurs for temperatures at 
which helium is neutral and thus makes little contribution to the 
opacity. The horizontal dotted lines indicate the mass limits of 
the models we calculated, so that the extent of the convective 
regions outside these bounds is unknown. These models are for a 
hydrogen envelope of mass Mjj = 10 — 6 Mq. 



Figure 10. Hydrogen Cooling Sequences 1: Luminosity 
Evolution: The solid lines indicate model white dwarf cooling 
sequences with a hydrogen envelope of 3 X 10~ 4 Mq for each of 
three representative total masses. The dashed lines are the equiv- 
alent sequences with a smaller hydrogen envelope of 10 _6 Mq. 
The difference in luminosities at earlier times is a result of the 
thicker hydrogen layer leading to a larger stellar radius (the ef- 
fective temperatures are closer - see Figure O) . 
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Figure 11. Hydrogen Cooling Sequences 2: Temperature 
Evolution: Here we show the effective temperature for the same 
sequences as in Figure lu. Note the large variation in temperature 
with envelope mass for the most massive models. This is the effect 
of the contribution of residual hydrogen burning at the base of 
the thicker hydrogen envelope. 



Figure 12. Helium Cooling Sequences: We show here pure 
helium models for the same masses as before. The evolution at 
effective temperatures below 6000 K is uncertain because of the 
inaccuracy of the photospheric opacities for neutral helium at 
these temperatures. 



Figure 13. V Band Cooling Sequences: We show here the 
absolute V magnitude determined from our cooling sequences. 
The curves are for 0.15, 0.25, 0.35 and 0.45 M Q . 



Figure 14. I Band Cooling Sequences: As for Figure |13[ but 
for absolute I magnitude. 



Figure 15. The Gravity- Effective Temperature Relation: 

The solid lines represent the true g — T e ft curves, and the dashed 
lines are the fits given by equations (|l])-(|])- 



Figure 16. The T = Mass-Radius Relation: The dotted 
line is the gravity as determined from the Paczynski (1967) mass- 
radius relation. The dashed line was obtained using the Hamada 
and Salpeter (1961) pure helium mass-radius relationjThe filled 
circles are for the thick H envelope models (equation | ) and the 
open circles, for the thin H envelope models (equation lOj). 
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